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Abstract 

The Feynman path integral representation of quantum theory is used 
in a non-parametric Bayesian approach to determine quantum potentials 
from measurements on a canonical ensemble. This representation allows 
to study explicitly the classical and semiclassical limits and provides a 
unified description in terms of functional integrals: the Feynman path 
integral for the statistical operator, and the integration over the space of 
potentials for calculating the predictive density. The latter is treated in 
maximum a posteriori approximation, and various approximation schemes 
for the former are developed and discussed. A simple numerical example 
shows the applicability of the method. 

PACS. 02.50.Tt Inference methods - 31.15.Kb Path-integral methods - 05.30.-d 
Quantum statistical mechanics 

1 Introduction 



The solution of the quantum many-body problem requires both techniques for 
solving the Schrodinger equation and knowledge of the underlying forces, often 
to be deduced from observational data. In the field of nuclear physics, forces are 
extracted from scattering data and ground state properties of the two~nucleon 
system, since no practicable basic theory of nuclear forces exists up to date. 
Given a data-based phenomenological nucleon-nucleon potential, one can, in 
principle, construct the related potential between two colliding nuclei. However, 
this is a formidable task which has been attacked only for a few simple cases in 
an approximate way, and a direct calculation of the nucleus-nucleus potential 
from observational data is highly desirable for practical applications like, e.g.. 
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in nuclear astrophysics. In solid state physics the basic force is known to be 
the Coulomb force, however, for a straight problem like the motion of a single 
electron under the influence of a crystal surface, one would prefer to deduce the 
respective potential directly from observational data rather than going through 
the full many-body problem of electrons and nuclei of the crystal. 

The reconstruction of such two-body forces or single-particle potentials from ex- 
perimental data constitutes a typical inverse problem of quantum theory. Such 
problems are notoriously ill-posed in the sense of Hadamard ^ and require ad- 
ditional a priori information to obtain a unique, stable solution. Well-knowm 
examples are inverse scattering theory and inverse spectral theory 3^ . They 
describe the kind of data which are necessary, in addition to a given spectrum, 
to identify the potential uniquely. For example, such data can be a second 
spectrum for different boundary conditions, knowledge of the potential on a 
half-interval or the phase shifts as a function of energy. However, neither a 
complete spectrum nor specific values of the potential or phase shifts at all en- 
ergies can be inferred by a finite number of measurements. Hence any practical 
algorithm for extracting two-body forces or single-particle potentials from ex- 
perimental data must rely on additional a priori assumptions like symmetries, 
smoothness, or asymptotic behaviour. If the available data refer to a system at 
finite temperature T ^ 0, one is led to the inverse problem of quantum statistics. 
In such a case, non-parametric Bayesian statistics 4 is especially well suited 
to include both observational data and a priori information in a flexible way. 

In a series of papers |5] , the Bayesian approach to inverse quantum statistics has 
been applied to reconstruct potentials (or two-body interactions) from particle- 
position measurements on a canonical ensemble. A priori information was im- 
posed through approximate symmetries (translational, periodic) or smoothness 
of the potential or by fixing the mean energy of the system. The likelihood 
model of quantum statistics (defining the probability for finding the particle at 
some position x for a system with potential V at temperature T) was treated 
in energy representation. 

In the present paper we apply the Feynman path integral representation of 
quantum mechanics jH] to calculate the statistical operator p and related quan- 
tities in coordinate space. This representation is of interest in the context of 
inverse problems in Bayesian statistics for two reasons: First, it allows to study 
the transition to the semiclassical and classical limits, relevant for example to 
atomic force microscopy so far treated on the level of classical mechanics. 
However, scales may soon be reached where the inclusion of quantum effects will 
be mandatory. Second, one obtains a unified description of Bayesian statistics 
in terms of path integrals. These are on one side the Feynman path integrals, 
needed in the likelihood model, and on the other side the functional integral 
over the space of potential functions V when calculating the predictive density 
as integral over the product of likelihood and posterior for all possible potentials. 

Our paper is organized as follows: An introduction to Bayesian statistics is 
presented in section 2, showing how Bayes' theorem about the decomposition 
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of joint probabilities can be used for the inverse problem of quantum statis- 
tics. A general expression is given for the likelihood of a quantum system with 
given potential V for a canonical ensemble, and the prior density is chosen as 
Gaussian process to implement a bias towards smoothness and/or periodicity 
of the potential V. This potential can be calculated from a non-linear differen- 
tial equation which results from the maximum posterior approximation for the 
predictive density. Two approximation schemes for solving the inverse problem 
of quantum statistics in path integral representation are developed. In the first 
variant (sections 3 and 4), the path integral in the likelihood is treated in sta- 
tionary phase approximation. The resulting stationarity equations with respect 
to the path are classical equations of motion for a particle in potential —V. 
These equations are to be solved simultaneously with the stationarity equations 
with respect to the potential, following from the maximum posterior approxi- 
mation in the above classical approximation of the path integral. In the second 
variant (section 5), the basic stationarity equations of the maximum posterior 
approximation are treated in terms of the Feynman integrals. These equations 
which involve the logarithmic derivatives of the statistical operator and the par- 
tition function are still exact in the sense of quantum theory. Approximation 
schemes for the statistical operator in coordinate representation and the corre- 
sponding partition function and their derivatives are developed in sections 6 and 
7. In section 6, the quadratic fluctuations for the statistical operator and the 
partition function around the classical paths of section 4 are determined, while 
section 7 deals with the respective derivatives. Three approximation schemes 
are studied under the general strategy that the statistical operator drops out in 
the logarithmic derivative. A simple numerical example is added in section 8 
to show that the path integral formalism can actually be used for problems of 
inverse quantum statistics. 



2 Bayesian Approach to Inverse Quantum 
Statistics 

The aim of this paper is to determine the dynamical laws of quantum sys- 
tems from measurements on a canonical ensemble. The method used is non- 
parametric Bayesian inference combined with the path integral representation 
of quantum theory which allows to study the transition to the classical limit. To 
be specific, we aim at reconstructing the potential V of the system from mea- 
surements of the position coordinate x of the particle for a canonical ensemble 
at temperature 1 / /3. 

The general Bayesian approach, tailored to the above problem, is based on two 
probability densities: 

1. a likelihood p{x\0, V) for the probability of outcome x when measuring 
observable O for given potential V, not directly observable, and 
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2. a prior density p {V) defined on the space V of possible potentials V. 



This prior gives the probabiUty for V before data have been collected. Hence it 
has to comprise all a priori information available for the potential, like symme- 
tries or smoothness. The need for a prior model, complementing the likelihood 
model, is characteristic for empirical learning problems which try to deduce a 
general law from observational data. 

These ingredients, likelihood and prior, are combined in Bayes theorem to define 
the posterior of V for given data D through 

P[xt\Ot) 

Eq. (2.1) is a direct consequence of the decomposition of the joint probability 
p {A, B) for two events A, B into conditional probabilities p {A\B) and p {B\A), 
respectively. Observational data D are assumed to consist of N pairs, 

D = {{xi, Oi)\l < i < N} , (2.2) 

where xt, Ot denote formal vectors with components Xi, Oi. Such data are 
also called training data, hence the label t- For independent data the likelihood 
factorizes as 

p{xt\Ot,V) = Y[p{xi\Oi,V) , (2.3) 

i 

where a chosen observable Oi may be measured repeatedly to give values Xi, 
equal or different among each other. The denominator in (2.1) can be viewed 
as normalization factor and can be calculated from likelihood and prior by in- 
tegration over V, 

p{xt\Ot) = J DVp{xt\Ot, V)p{V) . (2.4) 

The V-integral in eq. (2.4) stands for an integral over parameters, if we choose 
a parametrized space V of potentials, or for a functional integral over an infinite 
function space. 

To predict the results of future measurements on the basis of a data set D, one 
calculates according to the rules of probability theory the predictive density 

p{x\0, D) = j DV p{x\0, V)p{V\D) (2.5) 

which is the probability of finding value x when measuring observable O under 
the condition that data D are given. Here we have assumed that the probability 
of X is completely determined by giving potential V and observable O, and 
does not depend on training data Z?, p{x\0, V, D) = p(x\0, V), and that 
the probability for potential V given the training data D does not depend on 
observable O selected in the future, p{V\0, D) = p{V\D). 
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The integral (2.5) is high-dimensional in general and difficult to calculate in 
practice. Two approximations arc common in Baycsian statistics: The first 
one is an evaluation of the integral by Monte Carlo technique. The second 
one, which we will pursue in this paper, is the so called maximum a posteriori 
approximation. Assuming the posterior to be sufficiently peaked around its 
maximum at potential V* , the integral (2.5) is approximated by 

p{x\0,D) ^ p{x\0,V*) (2.6) 

where 

V* = argmaxygyp(F|i:>) = argmaxy g y p (xtIOt, V)p{V) (2.7) 

according to eq. (2.1) with the denominator independent of V. Maximizing 
the posterior p{V\D) with respect to ^ G V leads to solving the stationarity 
equations 

5vPiV\D) = = Sv{p{xt\Ot, V)piV)) (2.8) 

where Sy denotes the fimctional derivative S/SV. Equivalent to (2.8) and 
technically often more convenient is the condition for the log-posterior 

Sv In p {V\D) = = Sylnp {xt\Ot, V) + Svlrip {V) (2.9) 

which minimizes the energy E {V\D) = —lnp{V\D) and will be used in the 
following. 

A convenient choice for prior p {V) is a Gaussian process, 

p{V) ~ exp {-| {V - Vo\K\V - Vo)} = e-^^M (2.10) 

where 

r[v] = {V-Vo\K\V-Vo) = j dxdx' [v{x)-vo{x)] K {x, x') {v{x')-vo (x')], 

(2.11) 

assuming a local potential V {x, x') — v{x)6{x — x'). The mean Vq repre- 
sents a reference potential or template for V, and the real-symmetric, positive 
(semi-)dcfinite covariance operator (7 K)~^ acts on the potential, measuring the 
distance between V and Vq. The hyperparameter 7 is used to balance the prior 
against the likelihood term and is often treated in maximum a posteriori approx- 
imation or determined by cross-validation techniques. A bias towards smooth 
functions u {x) can be implemented hy K = —cP / dx^ choosing vq (x) = 0. If 
some approximate symmetry of v (x) is expected, like for a surface of a crystal 
deviating from exact periodicity due to point defects, one may implement a 
non-zero periodic reference potential vq {x) in eq. (2.11). 

The likelihood for our problem follows from the axioms of quantum theory: The 

probability to find value x when measuring observable O for a quantum system 
in a state described by a statistical operator p — p (V) is given by 

p{x\0,V) = Tv{Po{x)p{V)} (2.12) 
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where Pq {x) = S| \x, ^) (x, ^| projects on the space spanned by the orthonor- 
malizcd cigcnstates \x, of operator O with eigenvalue x, and the label ^ 
distinguishes degenerate eigenstates with respect to O. If the system is not 
prepared in an eigenstate of observable O, a quantum mechanical measurement 
will change the state of the system, i.e., will change p. Hence to perform re- 
peated measurements under same p requires the restoration of p before each 
measurement. For canonical ensembles at given temperature, 

_ eM-0H) 
^" Trexp(-/3i?) ^^-^^^ 

with Hamiltonian H = T -\- V and temperature 1 / /?, this means to wait 
between two consecutive observations until the system is thermalized again. 
Choosing the particle position operator x as observable O, the probability for 
value Xi is 

p{xi\x, v) = Z-'Tr{\xi) {Xi\e^p{-(3H)} = (2.14) 

with partition function 

Z =Trexp(-/3if) = j dx {x\exp {-p H)\x) (2.15) 

where we have dropped the label ^ to simplify notation. For N repeated mea- 
surements of f with results Xi,i = 1, ...N, one has under the above assumptions 
of independent measurements 

p{xt\Ot, V) = l[p{xi\x, v) [{xi\e-'^"\xi)Z-^] . (2.16) 

i i 

Combining eqs. (2.10), (2.11) and (2.16) leads to the posterior 

PiV\D) r. A. ^l[{x,\e-^"\x,)^ exp (-Jtm) = exp (-i? (y|D)) 

(2.17) 

with energy functional 

E{V\D) = -Y,^^{xi\e-^"\xi) + NlnZ + ^T[v] , (2.18) 

i 

functional F [v] defined in eq. (2.11). The corresponding stationarity equations 
(2.9) in explicit form 

f 4^(^de-^^N.) N 6Z J 5T 

f^^ {x,\Q-fi^\xi) Z 5v{x) 2 bv(x) ' ^ ' 

with 

^^=i^(.(x)-.o(x)), (2.20) 
determine the potential v [x). 
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In a series of papers 0, eqs. (2.19) have been studied successfully in energy 
representation for a variety of choices for prior p{V). This requires solving 
the Schrodinger equation H \(l)a) = Ea \(t>a), (</>a|</>/3) — ^ap, which allows to 
calculate the functional derivatives 6y Ea and Sy (j)a (cf. appendix) needed in eqs. 
(2.19). In the following sections we shall apply the path integral formulation 
of quantum mechanics in order to study the semiclassical as well as classical 
regimes. 



3 Likelihood in path integral representation 

The matrix elements appearing in eq. (2.16) can be written as path integrals 



j Dq (r) exp I 1 dr {dq / drf + v {q)) I . 



q(l3h)=Xi f [3h 

q(0}=Xi ^ 



(3.1) 

They are related to those of the time development operator of quantum mechan- 
ics by Wick rotation in the complex time plane. The corresponding variable 
transformation 

t = -ir (3.2) 
replaces real time t by imaginary time r and velocity dq / dt by 

dq .dq , . 

inducing a change of sign in the kinetic energy term. 

Representation (3.1) is understood as abbreviation of an infinite dimensional 
integral when dividing the interval [0, K\ into equidistant segments of length 
£ = (3Ti/ M, coding the path g (r) at discrete points — ekhy qk = q{Tk) 
and taking the limit M oo 

M ^ /M-l 

m \ "2" 



( M / 

e \ ^ / m 

q(ph)=x. 



h^\2 

fe=l \ 



, k=l / 

-v{qk)\ 



qk - Qk-i 



e 



with 



g(0)=Xi 

q{l3h) = Xi 



j Dq{T) exp , 



(3.4) 



m - 1 



9(0) = Xi 
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and Euclidean action 



/3?i 

m .2 



S[q]= I dT +v (q) (3.6) 







where we have introduced q = dq/dr for short. The boundary values are fixed, 
qo = qM = Xi. The partition function Z as trace in coordinate space can be 
written as path integral over all periodic functions q (r) of fixed period /3 h 

q{l3h)=x 

Z = Ti{e-f'") = j dx{x\e-^"\x) ^ j dx j (r) exp | - i [g] 



q(0)=x 



j 7^g(r) exp |-i5[(7]| 



(3.7) 



g(0) = ?(/3 h) 



The path integral for (x| exp (— /3i?)|a:), in lowest order stationary phase ap- 
proximation, is given by 

{x\e^^{-l3H)\x) =A., exp (^-^^S[q^]^ (3.8) 

with qx (t) being the solution of the classical equations of motion (|4.3|l with 
boundary conditions 

fe(0) = g.(/3?i) = a: , (3.9) 

and factor comprising the quadratic fiuctuations around the classical path 
qx (r) (see section 6, eqs. H6.7II . (j6.9|l V The x-integral in Z, eq. (3.7), can 
also be treated in stationary phase approximation. The action S [q] depends on 
X through the boundary values of q (r) , and its derivative with respect to the 

upper (lower) boundary value of q yields the corresponding momentum p. 

The stationarity condition for S thus poses the additional boundary condition 

p{(5h) - p{Q) ^ Q . (3.10) 

Hence one has to find xq such that the solutions of the classical equations of 
motion H4.3|l fulfill boundary conditions for both coordinate q (r) and velocity 

q{0) = q(/3h) = xq and q{0) = q{f3n). (3.11) 

Then 

Z = Aoexp (3.12) 

in lowest order stationary phase approximation, with Aq the analogue of A^, 
eq. ((231) • 
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4 Maximum posterior in stationary phase ap- 
proximation 



In the representation (|3.1|) . H3.7|l . the posterior density reads 

P{V\D) J (jj Dq, (r)j exp |-i F, fe, v] } (4.1) 



with total action 



fe, i;] = 5 te, i;] + ?i In Z M + 1^ r H , (4.2) 

inserting eqs. (3.4) and (3.6) into (2.17) and (2.18). Note that to each data 
point Xi is assigned its own path integral. 

Following the reasoning in section 2 for the w-integration we shall treat the 
integrals (3.4) in stationary phase approximation, looking for paths q (r) which 
minimize the action S [q] and account for the main contribution to the integrals. 
The corresponding stationarity equations, 

S S d 
Q = — = ~mq, + —V [qi) for i = I, ... , N (4.3) 
qi dq, 

with boundary conditions 

q,{0) ^ q,{ph) = X, (4.4) 

are the classical equations of motion for a fictitious particle of mass m in the in- 
verted potential —v (q) with boundary conditions determined by the data points 
Xi. Their solutions serve as starting points for a quantum mechanical expansion. 

For each path qi (r) the energy 

E^ = ^mq^ ~ v{q^) (4.5) 

is conserved. Equations H4.3() . (|4.4() have to be solved simultaneously with the 
stationarity equations (|2.19() . explicitly for F = F,;: 



5F 
Sv {x) 



J2 J drSiq.ir) - x) - phN{x\^^\x) 



+ -/h dx' K {x, x'){v{x') - vo{x')) 



(4.6) 



for the choice (2.10) of the prior p{V). For the derivative of In Z we refer to 
the appendix, eq. IjA.lOjl . 
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The integral in the first term of 1)4. over ^-distributions can be evaluated, for 
simple zeroes of the arguments, with the help of eq. (|4.5|l . 



n j! - J- 3i - ' 



n. ix) ^^^^^ 



rii (x) being the number of times tj. with qi (tj. ) = x, < Tj. < (3h. In the 
second term of (|4.5|l . the path integral for (a;| exp {—(] H)\x) is in lowest order 
stationary phase approximation given by (|3.8(l . and analogously Z by H3.12() . 

A compact and instructive form of condition 1)4. 6|l is obtained by multiplying 
with some arbitrary observable / (x) / N j3Ti and integrating over x, 

N . . ^-13 H 



= / drf{qdr)) - / dx f {x) {xf-^\x) 



+ ^ j dxdx'f (x) K{x, x') [v (x') - v„ (x')] (4.8) 

where /j denotes the mean of / with respect to (imaginary) time r along path 
qi (r) and (/) the thermal expectation value of observable /. Condition (|4.8() 
reminds of the ergodic theorem of statistical mechanics |H] concerning time and 
ensemble average, there are, however, differences in three respects: 



1. the time average in (|4.8|l is over a finite interval only, 

2. paths qi (r) refer to boundary conditions 1)4.4)1 rather than to initial con- 
ditions for q (r), q (r), and 

3. the prior gives a contribution to ()4.8)l . non-zero in general, in contradiction 
to the ergodic theorem. 



In the high temperature limit, (3^0, the prior term dominates condition 1)4. 8)1 . 
as expected, since the first two terms in ()4.8ll become /3-independent. Prior 
knowledge p {V) completely determines the maximum posterior solution. In 
contrast, the prior on v becomes negligeable at low temperature, corresponding 
to large /3- values, and the first two terms of ()4.8)l fulfill the ergodic theorem. In 
fact, the potential 

N 

V {x) = — lim a > S {x — Xi) (4-9) 

i = 1 
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is a solution of (|4.()|) . if the prior can be neglected. For the corresponding 

N 

classical potential —v (q) = + liina^oo i ^ (q ~ ^i) the equations of motion 

i = i 

(|4.3|) with boundary conditions (|4.4|l have unstable solutions 

q, (r) = a:, , (4.10) 

and the first term in 14.611 reads 

t"-' N N 



J2 J dT5{q,{T)-x) = / dT^5{x,-x) = (3h^5{x,-x) . (4.11) 

i— 1q pi— 1 2—1 

For large /3- values 

(x|e-^«|a:} ^ ^ (.tI^o^) (v'o,»e-'5^° (4^^2) 



where the A^-fold degenerate quantum ground state {x\ipoi) is strongly localized 
by potential v (x), eq. (|4.9|l . around the data points Xi such that 

\{x\ipo^)? ^ 5{x-x,) (4.13) 

N 

in proper normalization. Hence, with Z — q-PEq _ j^q-PEq^ second 

_ 4 = 1 

term in H4.5|l . 



i=l 



1 ^ 

Y.5{x-x,), (4.14) 



cancels the first one, eq. (j4.11|l . 



So far we have restricted ourselves to position measurements. If given data refer 
to other observables, one can use closure relations to calculate the required 
matrix elements in those observables while retaining the above path integral 
formalism. A typical example would be particle momenta rather than positions. 
In this case, Fourier transformation leads to 

{p\e-^"\p) = J dxdx' {p\x'){x'\e-l^"\x){x\p) 

^ J dxdx' J Dg(r) exp |-i (5[g] + ip(a;' - a;))| 

= J Dq{T) exp S^-^(^S[q] + ip{q{f3n) - qiO)))"^ (4.15) 

where the integration is over all path q (t) of the interval [0, ph]. Integral 
()4.15|l may then be calculated in saddle point approximation. Before presenting 
a numerical case study to show that the above formalism is feasible in practice, 
we will study an alternative approach to the path integral representation. 
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5 Maximum posterior stationarity equations in 
path integral representation 



In this section we start from the general stationarity equations (2.19) for the 
posterior, restricting ourselves to the case of position measurements for the sake 
of definiteness. The terms in (2.19) which result from the likelihood can all 
be expressed in terms of the basic path integral (3.1). As the matrix elements 
{xi\e~^ \xi) and partition function Z = J dx (a;|e~'^ need no further com- 
ment, we will directly proceed to evaluate their derivatives with respect to v. 

According to (3.7), the functional derivative of Z leads to 



SZ 



5v{x') h 

q{0)=q{l3h) 



I J Dq{T) exp j dT'5{q{T') - x') 



(3h 

IJdr' J Dq{T)cxp!^-^S[q]js{q{T') ^ x'), 



q{0)=q{l3h) 



(5.1) 



interchanging the order of integration in the second step. To evaluate the above 
path integral we observe that action S [q] and measure Dq (t) are invariant under 
cyclic shift of each path q (r) around some arbitrary value t'. This is displayed 
in Fig. 1 where the two paths qi (t), q2 (t) cover the same set of values in the 
interval [0, ph] and thus generate the same value for the integral J drv (q (r)). 
The same reasoning holds for the derivatives gi, q2, hence J dr {q (r))^ also has 
the same value for the two paths. Therefore the path integral in l|5.1|) . running 
over all periodic paths which go through point x' at time r', can be expressed 
as the path integral over all paths which start and end at x': 



q(0) = q{f3h) 
q(l3h) = x' 



Dq{T) exp |-i5[9]|(5(g(T') - x') 



j i?g(r) exp |-i^[g]| 



g(0) = x' 

= {x'\e-P"\x') , (5.2) 

using eq. (3.1). Note that (|5.2(l holds independent of the choice of r'. The 

pn 

remaining integral in (|5.1(l is then trivial, J dr' = f3h, confirming the expected 



result 

= r47T / dx{x\e-^"\x) = -P{x'\e-^"\x') . (5.3) 
V [x') V [x') J 
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In the functional derivative of matrix elements {xi\e ^^\xi) with respect to v, 

13% q(J3h)=Xi 

^-^(x,|e-^^|:r,) ^ -\ j dr' j Dq{r) l^-^-S [q]^ 5 {q{r') ~ x') 

q{0)=xi 
pn q(l3h) = Xi 

IJdr' I i?g(T) exp , (5.4) 



q{0)=x, 
q{r ) — X 



the path integral is split into two separate integrals according to Fig. 2. Taking 
into account the boundary conditions for the paths qi (r), (72 (''') on the r-axis, 
one obtains, under the r'-integral, a product of non-diagonal matrix elements 
of the statistical operator at different temperatures, 

5{x,\e-P"\x,) _ 
6 V (x') 

ph q2(ph)=Xi qi{T')=x' 

-IJdr' J Z?(72(r) expj-i^aNj J (r) exp j-^ 5i [gi]| 

92 (r') = x' qi (0) = Xi 

^ -\j rfT'(x,lexp |- (^/3 - i/||x')(x'|exp 




dp' (xi|exp{-(/3 - l3')H}\x') cxp {-/?' 

(5.5) 

where the full action S [q\ is split into parts, 

r' 

SM= jdT[^ql+v,{q,)) 


and (5.6) 



S2 [92] = j dr {^ql + V2 (<72 



The results (|5.3(l . H5.5|l and 1)5.6(1 obtained from the basic path integral (3.1) are 
still exact, in particular they are strictly equivalent to the respective expressions 
in energy representation (see appendix). With the above formulae, stationarity 
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equations (2.19) read 



N 

E 

i= 1 



(x,|exp(-(/3 - P')H)\x'){x'\expi-P'H)\x,) e'^" 

dp , , Nf3{x \^^\x 

{x.i\exp(-P H)\xt) Z 
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2Sv{x') 

for the prior of eqs. (2.10), (2.11). 



-jK{v{x') ~ vo{x'j) 

(5.7) 



In the foUowing sections we shaU study approximation schemes for the above 
derivatives H5.3() and 15. 5() . It wiU be shown under what assumptions the exact 
(quantum mechanical) resuh (|5.7(l approaches the approximate (semiclassical) 
form (|4.6|l . and how to find corrections to H4.6(l taking into account quantum 
fluctuations around the classical paths qi (r) of equations H4.8|l . H4.4|l . 



6 Quadratic fluctuations 

Matrix elements 

q{l3h)=x 

(xle-f^'^lx) ^ J i?g(T) exp (6.1) 

can be handled by standard techniques, starting from the stationarity equa- 
tions H4.3|l . H4.4|l . They read, dropping the label i and using v' = dv/ dqx for 
simplicity of notation, 

6 S 

= — = -mq + v' (q) with q(0) = q(j3h) = x . (6.2) 
Sq 

The stationary solutions qx (t) of 16.2|l yield the main contribution to the path 
integral (|6.1|) ; fluctuations around these solutions qx (r) can be taken care of by 
a variable transformation 

g(r) = q^ir) + r (t) with r (0) = r(/3n) = . (6.3) 

Assuming that only small deviations of q (r) from qx (t) are important for the 
integral ()6.1|l . we approximate 

v{q) = v{qx) + (q - Qx) v (qx) + ^ (? - Qxf v {Qx) (6.4) 
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and find for the action 

S[q] = J dr 


I3n 



/m .2 TO .2 . . / n\ 

(^y + y ^ + mq^ r + v{q,j, + r)j 



0n I3h 



I3h 



+ / dr {mq^r + v' (q^) rj 



(6.5) 



where the last term vanishes by virtue of H6.2|l . H6.3|l and partial integration 



dr [rnqxT + v' (q^) r] = [m q^ r 



(3h 



(6.6) 



For the additive action H6.5|) . H6.6|l the matrix element H6.1|l factorizes 

(a:|e-^^|x) - A,cxp|-i5[fe]| 



exp 



■ S [q^c] 



Dr (t) 



[3h 



exp { -- / dr -f^ + - w (q^) 



r (0) =0 

1 // 



(6.7) 



exp 



■ S [qx] 



Dr (t) 



r (0) =0 

S'S\q\ 



Sqir) 



•(r) 



with Hesse-matrix 
5^S[q] 



5q{T')5q{T) 



Hr - r') (-TO^ + v" {q,iT))j 



(6.8) 



in approximation 1)6. 4|l . For the integral in 1)6.7(1 to be well-defined, the Hesse- 
matrix 1)6. 8|l has to be positive-definite. Under the expansion ()6.4)l this holds if 
V (x) > for all x. 

The remaining path integral ()6.7)l can be evaluated by the 'shifting method'. 
We shall simply recall the result, known in the literature as van Vleck-formula 
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m 

(xle-'^^lx) = A,exp|-i5[g,]| 



f3h 

/ 

1 d^s[q,y ' 



exp (6.9) 



with (r) a solution of 

v" (q^ir)) = (6.10) 

Ka; (t) 

If (7a; does not vanish on the path (t), one can choose 

Kx (t) = qx (t) , (6.11) 

as is easily seen by differentiating H6.2|l with respect to t. Otherwise we look for 
a linear combination of the two linearly independent solutions k^x^ (t) — qx (t) 
and 

T 

The latter solution follows from the fact that the Wronski determinant of eq. 
(j6.2(l is constant. 

For the partition function Z we use the result H6.7II , (|6.9|l for the matrix element 
of the statistical operator 



Z = J dx' Ax' exp j--^ S* 



(6.13) 



and the x'-integration is done numerically. Combining H6.7|) and H6.13|l results 
in the normalized matrix element of the statistical operator in coordinate rep- 
resentation 

J dx' Ax' exp [-j- S [qx'l} 

For large masses m, formula l|6.7(l reproduces the result of classical statistical 
mechanics. In this case, the equations of motion (|6.2|l simplify, 

(jx ~ — v' {qx) — > for TO — > cxD , (6.15) 

TO 

and are solved by the static paths 

qx (r) = X (6.16) 
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for the boundary conditions of (|6.2I) . Then from H6.7() 

{x\e~'^"\x) ^ exp{-Pv{x)) J Dr {t) exp l-^ Jdr^^^r'^ + ^v" {x)r■ 
r(o)=o I 

(6.17) 

In the hmit of large masses, v (x) / m 0, the remaining integral in H6.17|l 
becomes independent of x so that the classical result is obtained: 



7 Matrix elements of the derivative of the sta- 
tistical operator with respect to the potential 

Three variants of approximations for the derivative 

ph q{ph) = Xi 



— 'xM'^" 



5 V {x') 



\x.,) = -\j dr' J Dq{T) exp {^^^m} 5{q{T')-x') 

q(0)=Xi 

(7.1) 

are presented. The general strategy is to find approximations such that in the 
logarithmic derivative of the statistical operator, needed in (2.19), the statistical 
operator drops out. 

In the first approach, we observe that the main contribution to the path 
integral stems from the stationary path q^^ (r), solution of eqs. (|4.^^() . H4.4|l . 
Hence the distribution 6 (q (t') — x') under the path integral in (|7.1(l may be 
replaced by 6 {q^^ (t') — x'), referring to the stationary path, in front of the 
path integral. In this approximation, 

fPh \q(pK) = Xi 

I dr'Siq^^ir') - x')\ J Dq (r) exp {-^^ m} , 



9(0)=x, 

the first term in the stationarity equation (2.19) takes the form 



(7.2) 



^^{x,\e-^"\x,) 1 7' 
{x^\e '^"\xi) h J 



in agreement with the first term in (|4.6|l . One may improve on the result ()7.3|l by 
applying approximation 1)7. 2|l after the classical action S [qx^] has been factored 
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out with the help of the variable transformation H().3|l . Then, using (|6.7|l . 

I3h r(f3h) = 



6 V (x') 



r(0)=0 



X ,5 (g,^ (r') + r (r') - x') exp <j 5 [r] 



/3?i r(/3?i) = 

"'"P^"^^''^"'^^ /'rfr'<5(g..(r') + r..(r')-a:') / Dr {r) e^p L^S [r 



r (0) =0 



fele-^^lx^ 



r(0)=0 

one finally has 

I3h 



n 

(7.4) 



where r^;. (r) is the solution of 

mr^i = w (gxi {T))rx, for r^, (0) = = r^;, {(3%) , (7.5) 

and 

-y^T Qmf2 + 1 (r))r2^ . (7.6) 



With lEHI), 

r(;3ft)=0 



)=exp|-i5fej| J i?r(r) exp |-i 5 h| , (7.7) 



A second possibility to evaluate the derivative of the statistical operator 
consists of inserting the spectral representation of the (5-distribution into H7.1|l . 
With g (t) = (t) + r (r) one obtains, in approximation (|6.4(l . for 

13 h 

j-l-^^{x,\e-^"\x.)= -1 exp |-i5[g,j|y'dry|^exp{zA(g,. {r') - x')} 



r([3h) = / f3h 

X y L)r (r)exp I -iydr + it)" (q^Jr^ - i A?i(5(t - r')r(T)^ 

r(0)=0 I 

(7.9) 

a path integral where a (5-like potential appears in the action in addition to the 
harmonic potential with time dependent frequency. Looking for saddle points 
of the path integral leads to the inhomogeneous equation of motion 

-m r.^^ r' A (t) + v"(fe, (r)) r^r, r' A = j A ?i (5 (t - r') (7.10) 
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with boundary conditions 



(7.11) 



Note that the solutions of (|7.10() . (|7.11(l are both A- and r'-dependent. The path 
fxiT'x is intimately related to the Green function R,^. (r, r') of the operator 
—1TI £^ + V {qxi (t)) with the above boundary conditions, 



-m i?,. (r, t') + v" (g,, (r)) i?,^ (r, r') = 5 (r - r') (7.12) 



with 



namely 



i?,, (0, t') = i?,. (/3?i, t') = 



(7.12' 



rx.r'A(r) ^ ihXRx^ (r, r') . (7.13) 

Furthermore, multiplying (|7.10() by r^^ t' \ (t) and integrating over r, one ob- 
tains for the quadratic part of the stationary action of H7.9|l 



l3h 



1 // 



iX , 



= — i?..(r',T')(7.14) 

Here use has been made of a partial integration, — Jdrr^ — j drrr, and of 
eqs. H7.in|l and H7.11|l . After further variable transformation. 



rir) = r^^r'X (t) + I (t) , 
one finds that in the action of eq. H7.9|l 



(7.15) 



/3h 



dT 



(m r^^ r' A (t) i (r) + v {q^^ (t)) r^. r' a (t) I (t) - ikXl (t) S (t ~ r')^ 



= [mr.,^r'x{r)l{T)fo'' = 
so that H7.9|l reads, under approximation (|7.2fl . 



(7.16) 



(5(rr,|e-^^|x,) ^ 1 



exp --^[g:^,] ^ /dT' 



/3fi 



dX 
27 



cxp< iA((7j;, (t') + 



i{0h) = o 



exp|-i5K.,,A]| y" DI{t) cxp S^-^S[l]^ 



(7.17) 
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with abbreviations S [vx. r' x] and 5* [I] defined according to 1)7. Using H7.14|l . 
the integral over A is of Gaussian type and can be carried out: 

^ exp \iX{qx. (t) + ^r^^r'xir') - x') 

^ exp L A (q^^ (r') -x')-^ R,, (r', r') \ (7.18) 



exp 



{q,^ (r') - xr 



V2 7r?ii?,. (r', t') I 2hR,^{T',T') 
The final result for the logarithmic derivative of {xi\e~^ ^ \xi) is 

' dr' ^ ^ , ^ , (7.19) 



(5 D (a;') 



{x,\e-f^"\xi) nj y/2nhR^^ (r', t' 

inserting 

(x,|e-'3^|x,) 



! (/3?i) = 

exp |-i5[9,.]|exp|-i^[r,.,,A]| ^ exp j-^ 5 [;] | (7.20) 

/(o) 

according to (|6.7(l . In comparison to 1)7. 3|l and 1)7.8(1 of the first approach, the 
(5-distributions in ((7.3(1 and ((7.8(1 are replaced in 1(7.19(1 by Gaussians which are 
normalized with respect to x' and whose widths are given by the Green functions 
R^. (t, r). Note that 

R.. {r\ r') =Jdr l^m(^^^^^Il^^ + \v" (,,^ [r)) R^ (r, A > 
^ ^ 

(7.21) 

for u" (q^.) > 0. Eq. ((7^ follows from ((71^ multiplied by R^, (r, r') and 
integrated over r. 

Finally, in our third approach we go back to eq. ((5.5(1 . The two matrix 
elements of the statistical operator under the /3'-integral can be expressed as 
path integrals separately. Stationarity of 

- / dp' {x,\exp{-{(3 - P') H)\x') {x'\exp{-(3' H)\x^) = 





/3 92 (/3 n) = Si qi{P'h)=x' 



jdp' j i?g2(T) exp |--i52[g2]| j i^gi (r) exp |-i 5i [gilj , 



q.2{P'h)=x' 5i(0) = Xi 



(7.22) 
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with respect to paths qi (r) and (72 (''") for [qi] and S2 [92] defined in H5.6|l . 
leads to the usual equations of motion and boundary conditions as indicated 
above. Stationarity with respect to (3' results in 

= ^ (^1 [qi] + S2 M) = jql iP' K)+v {q, [P' h))-^ qi (/?' n)-v {q^ (/?' Ti)) . 

(7.23) 

Since v (qi (/?' h)) — v (52 iP' ^)) by virtue of the boundary conditions, we have 
qi {/3' h) — ±q2 (/?' fi) for the velocities, and the energies Ej — ^niq^ ~ v (qj) 
are the same for both paths with j = I and j = 2. Our final result is 

{x,\e-0"\x.) = a:- [-% ['^^l + ['^^l ^ ^ t'^-l) 

(7.24) 

where qx^ is the solution of (|6.2|l with q^. (0) = Xi ^ qx^ (/3 K), Ap is a norming 
constant due to the stationary phase approximation of the /3'-integral in (|7.22() 
and v4i, ^2, stand for the fluctuations around the classical solutions as in 



8 Numerical case study 

In this section we present numerical results for a simple, one-dimensional model, 
which merely serve to demonstrate that the path integral technique can be used 
in actual practice within the Bayesian approach to inverse quantum statistics. 
We will discuss in turn the classical equations of motion (|4.3|l with boundary 
conditions l|4.4|l and the stationarity equations (|4.6(l of the maximum posterior 
approximation, which eventually have to be solved simultaneously. 

For a numerical implementation we discretize both the time t, parametrizing 
some classical path q (r), and the position coordinate x, upon which the poten- 
tial V {x) depends. The time interval [0, (3 fi] is divided into iir equal steps of 
length 

e = P jur (8.1) 

choosing units such that ?i = 1 . A path q (r) is then coded as vector q with 
components q^ = q{Tk) for Tk = e k; k = 0, 1, Ut- The potential v (x) 
is studied on an equidistant mesh of size Ux in space, choosing — tIt-. To 
match the equidistant values of coordinate x to the corresponding values of the 
classical path q (r) we may either round up or down the function values qu or 
linearly interpolate the potential between equidistant a:-values. 

In their discretized version, the classical equations of motion of our fictitious 
particle in potential —v (q) read 

7TL 

= -{qk + 1 - 2qk + qk~l) + v' (qk) ; fc = 1, 2, ... _ 1 , (8.2) 

£^ 

and are to be solved with boundary conditions 

qo = X = qn^ . (8.3) 
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Eqs. I|8.2|l and l|8.3|l amount to solving the matrix equation, for given v {q), 



= 



•■• 

1 -2 1 ••• 

1-210 



■•■ 
V ••• 

/ X \ 



91 



-2 1 



V / 



= Aq + t{q) 
which is done by iteration according to 



(8.4) 



V, {q^'^ + (■'"))) . (8.5) 



Step length ry^ in H8.5|l can be adapted during iteration. Having solved H8.4|l 
for various boundary values Xi, we can calculate the likelihood p {xi\Oi, V), eq. 
(2.14), in classical (eq. H6.18|l ') and semiclassical (eq. (|3.8|) ') approximations, 
and the exact quantum statistical result (see appendix). 



As example we consider a potential of the form 



V (x) 



1 



1 + cxp (i(|a; - 15| -4)) 



(8.6) 



on an equidistant mesh with = rir = 30, shown in Fig. 3., upper left part. 
The right hand side of Fig. 3 displays the potential — i; {q) together with the 
range and energy of solutions qx (t) of (|8.4II for various boundary values x (upper 
part), and the solutions qx (t) as functions of r. Note that solutions qx (t) refer 
to a boundary value problem in the fictitious potential —v (q) rather than to the 
initial value problem of classical mechanics in potential +v (q) . The probabilities 
p {xi\Oi, V), eq. (2.14), in the lower left part of Fig. 3 exhibit the difference of 
the classical and semiclassical approximations to the exact quantum statistical 
result. As expected on account of the uncertainty relation, the variance of the 
probability distribution increases when going from the classical limit to the exact 
quantum mechanical calculation. The 3 curves coincide in the classical result, 
if temperature or mass are increased. 

To evaluate the first term of stationarity equation (|4.6|l for one-dimensional 
models, one should not use eq. (|4.7|l : on every one-dimensional, periodic path 
the velocity takes the value zero for at least one value of r. A zero of the 
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argument of the (^-distribution at that value of t will not be a simple one. We 
have, therefore, for the discretized values Tfe, replaced the value q{Tk) by the 
nearest integer of the interval [0, n^], and the (5-distribution in 14.7(1 by the 
Kronecker symbol, hence 



drS [qi (r) - x) 



•7) 



The matrix elements of the second term {x\ exp (— /3 H)\x) are calculated semi- 
classically according to ((3.8(1 . In the prior (2.10) we use 



with 



/ -2 1 





-2 1 



\ 



V 



1 -2 
1 



(8.9) 



2/ 



thus demanding smoothness for the potential v {x) to be determined. In our 
actual calculation we have sampled = 15 data from the discretized version 
{rix = 30) of potential 



V (x) 







1 for a; e [5, 25] 



elsewhere 



(8.10) 



for P = 10. Eq. 14.6(1 is then solved, simultaneously with eqs. 18.4(1 . 18.5(1 . by 
iteration, using the gradient descent algorithm. The hyperparameter 7 is chosen 
such that the depths of the reconstructed potential and the true potential ((8.10(1 
are approximately equal as shown on the left hand side of Fig. 4. The right 
hand side shows the empirical density of data together with the likelihood for 
the true potential and for the classical, semiclassical and quantum mechanical 
reconstructions. For sufhciently heavy masses the gross shape of the potential 
is recognized; classical, semiclassical and quantum mechanical likelihoods are 
approximately the same. With decreasing mass, the differences of classical, 
semiclassical and quantum mechanical likelihoods become more pronounced, 
with the double-hump structure of the potential still recognized (Fig. 5). To 
better reproduce the absolute value of the potential minima one may decrease 
parameter 7 at the expense of distorting the symmetrical shape of the potential, 
like in Fig. 4. 
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9 Conclusion 



In this paper wc have developed the inverse problem of quantum statistics in 
path integral representation which supplements the energy representation used 
in a number of recent publications. The advantage of the path integral rep- 
resentation in this context turns out to be twofold: First, one can study the 
semiclassical and classical limits which are of interest for the analysis of exper- 
imental data as obtained from atomic force microscopy. Second, with the path 
integral representation for the likelihood and the functional integration over pos- 
sible potential fields, one obtains a unified description for the basic equations 
of Bayesian inverse quantum statistics. Various approximation schemes have 
been studied for calculating, in this representation, the statistical operator and 
its derivatives which are the essential quantities in maximum posterior approx- 
imation. In particular, the classical limit is obtained and quadratic quantum 
fluctuations are calculated. A simple numerical example is presented to demon- 
strate the actual applicability of this approach which is expected to be useful for 
analysing experimental data when spatial distances are reached which resolve 
nanostructures like in recent atomic force microscopy. 



Appendix 

Matrix elements (.T'|e~''^|3;} of the statistical operator and their fimctional 
derivatives with respect to V, needed in eq. (2.19), are easily calculated in 
energy representation. We start from the Schrodinger equation 

H \(bo.) ^ [T + V) = \(^^) , (A.l) 

together with orthonormality and closure of eigenfunctions, 

^ |<^„)(.^„| = n. (A.2) 
a 

For the derivatives of 

{x'\e-^"\x) = ^{x'\<l>c.)e-^^''{<p^\x) 

a 

= Y.^M^-{x')e-^^" (A.3) 

a 

and of 

Z = j dx{x\e-l^"\x) = I dx^{x\(l>c.)e-^^- {(l)a\x) 

■' ■'a 

= 5] / rfa; {x) {x) e-" ^° = ^ e-^ ^" (A.4) 
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we need d Ea / Sv (x) and S (pa {x') / 6v (x) . These derivatives are obtained by 
variation of the Schrodinger equation, 

^ + H \Sv = ^ + Ea \5v M (A.5) 
V V 

where, for a local potential V {x' , x ) — v (x ) S {x' ~ x ), 
, SH _ 6H{x', x") 



(X - — X } = 



Sv{x) Sv{x) 

in short 



^ S{x' - X )S{x - X ) , (A.6) 



= 1^) {x\ . (A.7) 

V [x) 



Multiplying (A.5) by {(j>a\ from the left and using the adjoint of (A.l), 

{^a\H = {^a\Ea , (A.8) 

yields with normalization (A. 2) 

ov (X) oV 

= {<j)a\x) {X\(ba) = \^a.ix)\^ ■ (A.9) 

Hence our first result, in agreement with H5.^{(l . is 

= -^^e-''^"/^ = -/?^e-^^=|<^„(x)p 
ovix) ^-^ ovix) ^-^ 

= -/3 ^(a;|0„)e-^^° = -/3 (xje"'' , (A.IO) 

a 

with closure (A. 2). 

To find the derivative of (pa {x') with respect to v (x) we rewrite (A.5) as inho- 
mogeneous equation, 

{Ea H) \5y Cba) = " ^) ■ (A.ll) 

Obviously all orbitals \(j)a') with Ea' = Ea are in the null space of the operator 
{Ea — H), hence {Ea ~ H) is not invertible in full Hilbert space. However, the 
Moore-Penrose method of the pseudo- inverse can be applied to solve (A.ll) for 
|5y 0q). The solvability condition states that the right hand side of (A.ll) has 
no component in the null space of {Ea — H). This is easily verified, multiplying 
(A.ll) by {(Pa' \ with£;„. = Ea-. 

{(Pa'\Ea - H\Sv^a) = = {(f>a' \^ - ^\cpa) , (A.12) 

V V 

applying H to the left according to (A.8). It is easy to control that 



^7 
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is the pseudo- inverse of {£„ — H), fulfilling the condition 

G„ - ) G„ = G„ . (A. 14) 

To obtain a unique solution of (A. 11), or equivalcntly (A. 5), we demand that 
\6v (pa) has no component in the null space of (Ea — H), 

{(t>a'\Sv (pa) = for all (pa' with Ea' = Ea ■ (A.15) 

This corresponds to fixing norm and phase of the eigenstates \(pa) and, in case of 
degeneracy, uses the freedom to work with arbitrary orthonormal linear combi- 
nations of the respective eigenstates. Applying the pseudo-inverse Ga to (A. 11) 
and using the orthonormality (A. 2) we thus find in the subspace where {Ea — H) 
is invertible 

\5v(l>a)= {Ea - E^)-^\cP^){<p^\5vH\(l>a) , (A.16) 

7 

or explicitly in coordinate representation with eq. (A. 6) 

7 

-E-y / -Eq 

With the derivatives (A. 9) and (A. 17), given in terms of the solutions of 
the Schrodinger equation, we can now construct the functional derivative of 
(a;'|e''^|a;) with respect tov{x ) according to (A. 3): 

— ^ {x'\e-^"\x) = -/3e-^^° ^ C(^) M^') \<Pa {x')? 
dv(x I '-^ 

\ / a 

(A.18) 



+ E {'^"^"'^ '^^^'"'^ "^^^^"^ '^"^^"^ + '^^^'"^ '^^^^"^ (P*a{x")<Pa{x')} ■ 



a, 7 

" #7 



In particular, for diagonal elements with x' = x 

S 



^^^-{x\e-^"\x) = -fie-P'^- \<Pa{x)?\(Pa{x")\ 

^ cxp(-/3-B„) / N , / N , t "\\ (^-19) 

+ E E - E l^" ^ '^'^ 



a, 7 
ai^l 
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This result, eq. (A. 19), is identical with (|5.5|l as can be shown by carrying 
out the /^'-integration in H5.5() using the energy representation of the matrix 
elements under the integral: 



p 

dv(x ) J 



(A.20) 



P for a = 7 , or else 



inserting the closure relation (A. 2). With the /^'-integration carried out, 
n' i-R' ij^. - w - ) 

1 _ exp(-/3(£;^ - E^)) 

y Ej — Eq, Ej — Ea 

(A.21) 

one obtains 

Sixle-f^^lx) R^-PE^STij, I \|2u I 

- -fie ' " |0a (a;)| |0a (a; ) 



Sv {x") 



a, 7 



q. e. d. 
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Figure 1: Equivalent paths qi (r) and 52 (t) in the interval [0, hf3], explaining 
eq. H5.2|l . (In the figure r is denoted t.) 
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Figure 2: Splitting path 9 (r) into parts qiir) and 92 ('''), explaining the 
transition from eq. H5.4|l to (|5.5(l . (In the figure r, r' are denoted t, t' , re- 
spectively.) 
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Figure 3: Comparison of classical, semiclassical and quantum mechanical like- 
lihood. Upper left part: original potential v{x). Upper right part: potential 
u (x) — —V (x), to be used in (|8.2|l : thin horizontal lines indicate range and en- 
ergy of paths Qx (t). Lower left part: classical (dotted line), semiclassical (full 
line) and quantum mechanical (dashed line) likelihoods. Lower right part: paths 
Qx (t) for various x-values. Parameters used are m = 0.1, f3 — 6, rir — 30, 
Ux — 30. (In the figure r is denoted t.) 
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Figure 4: Bayesian reconstruction of potentials using the path integral method. 
Left part: original potential (thin line) and reconstructed potential (thick line). 
Right part: relative frequencies of sampled data (bars), likelihood of the true 
(thin line) potential and of the reconstructed potential: semiclassical density ob- 
tained by iteration (thick line), classical (dotted line) and quantum mechanical 
density (dashed line). Parameters: /3 = 10, m = 1, 7 = 5, TV = 15. 
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Figure 5: Bayesian reconstruction of potentials, using the path integral method, 
for small masses: Graphics as in Fig. 4. Parameters: /3 = 10, m = 0.05, 
7 = 10, iV = 15. 
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